1 Abstract

The muHVT package is a collection of R functions to facilitate building topology preserving maps for rich multivariate data analysis. Tending towards a big data preponderance, a large number of rows. A collection of R functions for this typical workflow is organized below:

1. Data Compression: Vector quantization (VQ), HVQ (hierarchical vector quantization) using means or medians. This step compresses the rows (long data frame) using a compression objective.

2. Data Projection: Dimension projection of the compressed cells to 1D,2D or 3D with the Sammons Non-linear Algorithm. This step creates topology preserving map (also called an embedding) coordinates into the desired output dimension.

3. Tessellation: Create cells required for object visualization using the Voronoi Tessellation method, package includes heatmap plots for hierarchical Voronoi tessellations (HVT). This step enables data insights, visualization, and interaction with the topology preserving map useful for semi-supervised tasks.

4. Prediction: Scoring new data sets and recording their assignment using the map objects from the above steps, in a sequence of maps if required.

2 Data Compression

Compression is a technique used to reduce the data size while preserving its essential information, allowing for efficient storage and decompression to reconstruct the original data. While Vector quantization (VQ) is a technique used in data compression to represent a set of data points with a smaller number of representative vectors. It achieves compression by exploiting redundancies or patterns in the data and replacing similar data points with representative vectors.

This package offers several advantages for performing data compression as it is designed to handle high-dimensional data more efficiently. It provides a hierarchical compression approach, allowing multi-resolution representation of the data. The hierarchical structure enables efficient compression and storage of the data while preserving different levels of detail. HVT aims to preserve the topological structure of the data during compression.Spatial data with irregular shapes and complex structures in high-dimensional data can contain valuable information about relationships and patterns. HVT seeks to capture and retain these topological characteristics, enabling meaningful analysis and visualization.This package employs tessellation to divide the compressed data space into distinct cells or regions while preserving the topology of the original data. This means that the relationships and connectivity between data points are maintained in the compressed representation.

This package can perform vector quantization using the following algorithms-

• Hierarchical Vector Quantization using \(k-means\)
• Hierarchical Vector Quantization using \(k-medoids\)

2.1 Hierarchical Vector Quantization

2.1.1 Using k-means

1. The k-means algorithm randomly selects k data points as initial means.
2. k clusters are formed by assigning each data point to its closest cluster mean using the Euclidean distance.
3. Virtual means for each cluster are calculated by using all datapoints contained in a cluster.

The second and third steps are iterated until a predefined number of iterations is reached or the clusters converge. The runtime for the algorithm is O(n).

2.1.2 Using k-medoids

1. The k-medoids algorithm randomly selects k data points as initial means out of the n data points as the medoids.
2. k clusters are formed by assigning each data point to its closest medoid by using any common distance metric methods.
3. Virtual means for each cluster are calculated by using all datapoints contained in a cluster.

The second and third steps are iterated until a predefined number of iterations is reached or the clusters converge. The runtime for the algorithm is O(k * (n-k)^2).

These algorithm divides the dataset recursively into cells using \(k-means\) or \(k-medoids\) algorithm. The maximum number of subsets are decided by setting \(n_cells\) to, say five, in order to divide the dataset into maximum of five subsets. These five subsets are further divided into five subsets(or less), resulting in a total of twenty five (5*5) subsets. The recursion terminates when the cells either contain less than three data point or a stop criterion is reached. In this case, the stop criterion is set to when the cell error exceeds the quantization threshold.

The steps for this method are as follows:

1. Select k(number of cells), depth and quantization error threshold.
2. Perform quantization (using \(k-means\) or \(k-medoids\)) on the input dataset.
3. Calculate quantization error for each of the k cells.
4. Compare the quantization error for each cell to quantization error threshold.
5. Repeat steps 2 to 4 for each of the k cells whose quantization error is above threshold until stop criterion is reached.

The stop criterion is when the quantization error of a cell satisfies one of the below conditions:

• reaches below quantization error threshold.
• there are less than three data points in the cell.
• the user specified depth has been attained.

The quantization error for a cell is defined as follows:

\[QE = \max_i(||A-F_i||_{p})\]

where

• \(A\) is the centroid of the cell
• \(F_i\) represents a data point in the cell
• \(m\) is the number of points in the cell
• \(p\) is the \(p\)-norm metric. Here \(p\) = 1 represents L1 Norm and \(p\) = 2 represents L2 Norm

2.1.3 Quantization Error

Let us try to understand quantization error with an example.

Figure 1: The Voronoi tessellation for level 1 shown for the 5 cells with the points overlayed

An example of a 2 dimensional VQ is shown above.

In the above image, we can see 5 cells with each cell containing a certain number of points. The centroid for each cell is shown in blue. These centroids are also known as codewords since they represent all the points in that cell. The set of all codewords is called a codebook.

Now we want to calculate quantization error for each cell. For the sake of simplicity, let’s consider only one cell having centroid A and m data points \(F_i\) for calculating quantization error.

For each point, we calculate the distance between the point and the centroid.

\[ d = ||A - F_i||_{p} \]

In the above equation, p = 1 means L1_Norm distance whereas p = 2 means L2_Norm distance. In the package, the L1_Norm distance is chosen by default. The user can pass either L1_Norm, L2_Norm or a custom function to calculate the distance between two points in n dimensions.

\[QE = \max_i(||A-F_i||_{p})\]

Now, we take the maximum calculated distance of all m points. This gives us the furthest distance of a point in the cell from the centroid, which we refer to as Quantization Error. If the Quantization Error is higher than the given threshold, the centroid/ codevector is not a good representation for the points in the cell. Now we can perform further Vector Quantization on these points and repeat the above steps.

Please note that the user can select mean, max or any custom function to calculate the Quantization Error. The custom function takes a vector of m value (where each value is a distance between point in n dimensions and centroids) and returns a single value which is the Quantization Error for the cell.

If we select mean as the error metric, the above Quantization Error equation will look like this:

\[QE = \frac{1}{m}\sum_{i=1}^m||A-F_i||_{p}\]

3 Data Projection

Projection mainly involves converting data from its original form to a different space or coordinate system while preserving certain properties of it. By projecting data into a common coordinate system, spatial relationships, distances, areas, and other spatial attributes can be accurately measured and compared.

HVT performs projection as part of its workflow to visualize and explore high-dimensional data. The projection step in HVT involves mapping the compressed data, represented by the hierarchical structure of cells, onto a lower-dimensional space for visualization purposes, as human perception is more suited to interpreting information in lower-dimensional spaces.Users can zoom in/out, rotate, and explore different regions of the projected space to gain insights and understand the data from different perspectives.

Sammon’s projection is an algorithm used in this package to map a high-dimensional space to a space of lower dimensionality while attempting to preserve the structure of inter-point distances in the projection. It is particularly suited for use in exploratory data analysis and is usually considered a non-linear approach since the mapping cannot be represented as a linear combination of the original variables. The centroids are plotted in 2D after performing Sammon’s projection at every level of the tessellation.

Denoting the distance between \(i^{th}\) and \(j^{th}\) objects in the original space by \(d_{ij}^*\), and the distance between their projections by \(d_{ij}\). Sammon’s mapping aims to minimize the below error function, which is often referred to as Sammon’s stress or Sammon’s error.

\[E=\frac{1}{\sum_{i<j} d_{ij}^*}\sum_{i<j}\frac{(d_{ij}^*-d_{ij})^2}{d_{ij}^*}\]

The minimization of this can be performed either by gradient descent, as proposed initially, or by other means, usually involving iterative methods. The number of iterations need to be experimentally determined and convergent solutions are not always guaranteed. Many implementations prefer to use the first Principal Components as a starting configuration.

4 Tessellation

A Voronoi diagram is a way of dividing space into a number of regions. A set of points (called seeds, sites, or generators) is specified beforehand and for each seed, there will be a corresponding region consisting of all points within proximity of that seed. These regions are called Voronoi cells. It is complementary to Delaunay triangulation is a geometrical algorithm used to create a triangulated mesh from a set of points in a plane which has the property that no data point lies within the circumcircle of any triangle in the triangulation. This property guarantees that the resulting cells in the tessellation do not overlap with each other.

By using Delaunay triangulation, HVT can achieve a partitioning of the data space into distinct and non-overlapping regions, which is crucial for accurately representing and analyzing the compressed data.Additionally, the use of Delaunay triangulation for tessellation ensures that the resulting cells have well-defined shapes, typically triangles in two dimensions or tetrahedra in three dimensions.

The hierarchical structure resulting from tessellation preserves the inherent structure and relationships within the data. It captures clusters, subclusters, and other patterns in the data, allowing for a more organized and interpretable representation. The hierarchical structure reduces redundancy and enables more compact representations.

Tessellate: Constructing Voronoi Tesselation

In this package, we use sammons from the package MASS to project higher dimensional data to a 2D space. The function hvq called from the HVT function returns hierarchical quantized data which will be the input for construction of the tessellations. The data is then represented in 2D coordinates and the tessellations are plotted using these coordinates as centroids. We use the package deldir for this purpose. The deldir package computes the Delaunay triangulation (and hence the Dirichlet or Voronoi tessellation) of a planar point set according to the second (iterative) algorithm of Lee and Schacter. For subsequent levels, transformation is performed on the 2D coordinates to get all the points within its parent tile. Tessellations are plotted using these transformed points as centroids. The lines in the tessellations are chopped in places so that they do not protrude outside the parent polygon. This is done for all the subsequent levels.

5 Prediction

Prediction basically refers to the process of making predictions or estimating future values or outcomes based on existing data patterns.In data prediction, a model is developed based on historical data or a training dataset, and this model is then used to make predictions on new, unseen data. The model captures the underlying patterns, trends, and relationships present in the training data, allowing it to make informed predictions on similar or related data points.

In this package, we use predictHVT function to predict each point in the test dataset.

Prediction Algorithm

The prediction algorithm recursively calculates the distance between each point in the test dataset and the cell centroids for each level. The following steps explain the prediction method for a single point in the test dataset:

1. Calculate the distance between the point and the centroid of all the cells in the first level.
2. Find the cell whose centroid has minimum distance to the point.
3. Check if the cell drills down further to form more cells.
4. If it doesn’t, return the path. Or else repeat steps 1 to 4 till we reach a level at which the cell doesn’t drill down further.

6 Example I: muHVT with the Torus dataset

In this section, we will see how we can use the package to visualize multidimensional data by projecting them to two dimensions using Sammon’s projection and further used for scoring

Data Understanding

First of all, let us see how to generate data for torus. We are using a library geozoo for this purpose. Geo Zoo (stands for Geometric Zoo) is a compilation of geometric objects ranging from three to 10 dimensions. Geo Zoo contains regular or well-known objects, eg cube and sphere, and some abstract objects, e.g. Boy’s surface, Torus and Hyper-Torus.

Here, we will generate a 3D torus (a torus is a surface of revolution generated by revolving a circle in three-dimensional space one full revolution about an axis that is coplanar with the circle) with 9000 points.

Raw Torus Dataset

The torus dataset includes the following columns:

• x: This column represents the X-coordinate of each point in the torus.
• y: This column represents the Y-coordinate of each point in the torus.
• z: This column represents the Z-coordinate of each point in the torus.

Lets, explore the raw torus dataset containing 12000 points. For the sake of brevity we are displaying first 6 rows.

set.seed(240)
# Here p represents dimension of object
# n represents number of points
torus <- geozoo::torus(p = 3,n = 12000)
torus_df <- data.frame(torus$points)
colnames(torus_df) <- c("x","y","z")

torus_df1 <- torus_df %>% round(4)
colnames(torus_df1) <- c("x","y","z")
torus_df1$Row.No <- as.numeric(row.names(torus_df))
torus_df1 <- torus_df1 %>% dplyr::select(Row.No,x,y,z)
Table(head(torus_df1))

We will first split the torus data into train and test. We will randomly select 9000 data points as training and remaining 3000 data points as testing data.

set.seed(42)
train_indices <- sample(1:nrow(torus_df), 9000)
trainTorus <- torus_df[train_indices, ]
trainTorus_data <- trainTorus %>% round(4)
test_indices <- setdiff(1:nrow(torus_df), train_indices)
testTorus <- torus_df[test_indices, ]

Raw Training Dataset

First of all, we will see the randomly selected training data containing (9000 data points). For the sake of brevity we are displaying first six rows.

trainTorus_data$Row.No <- as.numeric(row.names(trainTorus_data))
trainTorus_data <- trainTorus_data %>% dplyr::select(Row.No,x,y,z)
row.names(trainTorus_data) <- NULL
Table(head(trainTorus_data))

Now let’s have a look at structure and summary of the training data.

str(trainTorus_data)
#> 'data.frame':    9000 obs. of  4 variables:
#>  $ Row.No: num  10801 2369 5273 9290 1252 ...
#>  $ x     : num  -0.686 0.047 1.415 0.245 -2.087 ...
#>  $ y     : num  -0.8709 -1.4714 0.0936 1.1402 0.0771 ...
#>  $ z     : num  0.454 0.849 0.814 -0.552 0.996 ...

summary(trainTorus_data)
#>      Row.No            x                   y                   z            
#>  Min.   :    1   Min.   :-2.997700   Min.   :-2.995600   Min.   :-1.000000  
#>  1st Qu.: 2988   1st Qu.:-1.151025   1st Qu.:-1.118100   1st Qu.:-0.716225  
#>  Median : 5986   Median : 0.022200   Median :-0.000600   Median : 0.016950  
#>  Mean   : 5988   Mean   :-0.002215   Mean   : 0.002805   Mean   : 0.004401  
#>  3rd Qu.: 8974   3rd Qu.: 1.140325   3rd Qu.: 1.125900   3rd Qu.: 0.719875  
#>  Max.   :12000   Max.   : 2.998100   Max.   : 2.999300   Max.   : 1.000000

Raw Testing Dataset

Now, lets have a look at randomly selected testing dataset containing(3000 data points).For the sake of brevity we are displaying first six rows.

test_dataset <- testTorus
test_dataset1 <- round(test_dataset,4)
test_dataset1$Row.No <- row.names(test_dataset)
test_dataset1 <- test_dataset1 %>% dplyr::select(Row.No,x,y,z) 
rownames(test_dataset1) <- NULL
Table(head(test_dataset1))

Now let’s have a look at structure and summary of the test data.

str(test_dataset)
#> 'data.frame':    3000 obs. of  3 variables:
#>  $ x: num  -1.482 0.792 -2.379 -0.843 -2.649 ...
#>  $ y: num  0.923 -1.348 1.799 -0.544 -0.574 ...
#>  $ z: num  0.9672 -0.8998 -0.1878 0.0755 0.704 ...

summary(test_dataset)
#>        x                    y                  z            
#>  Min.   :-2.9976672   Min.   :-2.99934   Min.   :-1.000000  
#>  1st Qu.:-1.1408711   1st Qu.:-1.09877   1st Qu.:-0.700378  
#>  Median :-0.0670732   Median : 0.06562   Median : 0.012098  
#>  Mean   : 0.0008702   Mean   : 0.03297   Mean   : 0.004486  
#>  3rd Qu.: 1.1404037   3rd Qu.: 1.14810   3rd Qu.: 0.713435  
#>  Max.   : 2.9995467   Max.   : 2.98818   Max.   : 0.999999

Now let’s try to visualize the torus (donut) in 3D Space.

knitr::include_graphics('torus_donut.png')

Figure 2: 3D Torus

Note: The steps of compression, projection, and tessellation are iteratively performed until a minimum compression rate of 80% is achieved. Once the desired compression is attained, the resulting model object is used for scoring using the predictHVT() function

In this section all the outlined workflow steps provided in the abstract section (Compression, Projection, Tessellation and Prediction) are executed at level 1.

HVT(
  dataset,
  min_compression_perc,
  n_cells,
  depth,
  quant.err,
  projection.scale,
  normalize = T,
  distance_metric = c("L1_Norm", "L2_Norm"),
  error_metric = c("mean", "max"),
  quant_method = c("kmeans", "kmedoids"),
  diagnose = TRUE,
  hvt_validation = FALSE,
  train_validation_split_ratio = 0.8
)

set.seed(240)
hvt.torus <- muHVT::HVT(
  torus_df,
  n_cells = 100,
  depth = 1,
  quant.err = 0.1,
  projection.scale = 10,
  normalize = F,
  distance_metric = "L1_Norm",
  error_metric = "max",
  quant_method = "kmeans"
)

compressionSummaryTable(hvt.torus[[3]]$compression_summary)

set.seed(240)
hvt.torus2 <- muHVT::HVT(
  torus_df,
  n_cells = 300,
  depth = 1,
  quant.err = 0.1,
  projection.scale = 10,
  normalize = F,
  distance_metric = "L1_Norm",
  error_metric = "max",
  quant_method = "kmeans"
)

compressionSummaryTable(hvt.torus2[[3]]$compression_summary)

set.seed(240)
hvt.torus3 <- muHVT::HVT(
  torus_df,
  n_cells = 900,
  depth = 1,
  quant.err = 0.1,
  projection.scale = 10,
  normalize = F,
  distance_metric = "L1_Norm",
  error_metric = "max",
  quant_method = "kmeans"
)

compressionSummaryTable(hvt.torus3[[3]]$compression_summary)

6.2 Step 2: Data Projection

The function sammonsProjection() utilizes the sammons function from the MASS package being called in HVT. Sammon’s projection is an algorithm that maps a high-dimensional space to a space of lower dimensionality while attempting to preserve the structure of inter-point distances in the projection.The centroids are plotted in 2D after performing Sammon’s projection at every level of the tessellation.

6.2.1 Iteration 1:

lets view the projected 2D coordinates after performing sammon’s projection on the compressed data for the first iteration where we set n_cells parameter as 100. For the sake of brevity we are displaying first six rows.

hvt_torus_coordinates <-hvt.torus[[2]][[1]][["1"]]
centroids <<- list()
  coordinates_value <- lapply(1:length(hvt_torus_coordinates), function(x){
    centroids <-hvt_torus_coordinates[[x]]
    coordinates <- centroids$pt
  })
centroid_coordinates<<- do.call(rbind.data.frame, coordinates_value)  
colnames(centroid_coordinates) <- c("x_coord","y_coord")
centroid_coordinates$Row.No <- as.numeric(row.names(centroid_coordinates)) 
centroid_coordinates <- centroid_coordinates %>% dplyr::select(Row.No,x_coord,y_coord)
centroid_coordinates1 <- centroid_coordinates %>% data.frame() %>% round(4)
Table(head(centroid_coordinates1), scroll = T, limit = 20)

Row.No	x_coord	y_coord
1	15.4686	9.1562
2	-12.3060	-3.5491
3	-6.9791	19.6759
4	9.5694	-0.5423
5	24.8946	17.7822
6	24.0559	6.8543

Lets see the projected Sammons 2D onto a plane with n_cell set to 100 in first iteration.

ggplot(centroid_coordinates1, aes(x_coord, y_coord)) +
  geom_point(color = "blue") +
  labs(x = "X", y = "Y")

Figure 3: Sammons 2D Plot for 100 cells

6.2.2 Iteration 2:

lets view the projected 2D coordinates after performing sammon’s projection on the compressed data for the Second iteration where we set n_cells parameter as 300. For the sake of brevity we are displaying first six rows.

hvt_torus_coordinates <-hvt.torus2[[2]][[1]][["1"]]
centroids <<- list()
  coordinates_value <- lapply(1:length(hvt_torus_coordinates), function(x){
    centroids <-hvt_torus_coordinates[[x]]
    coordinates <- centroids$pt
  })
centroid_coordinates<<- do.call(rbind.data.frame, coordinates_value)  
colnames(centroid_coordinates) <- c("x_coord","y_coord")
centroid_coordinates$Row.No <- as.numeric(row.names(centroid_coordinates)) 
centroid_coordinates <- centroid_coordinates %>% dplyr::select(Row.No,x_coord,y_coord)
centroid_coordinates2 <- centroid_coordinates %>% data.frame() %>% round(4)
Table(head(centroid_coordinates2), scroll = T, limit = 20)

Row.No	x_coord	y_coord
1	23.7284	5.0557
2	-11.2747	1.3672
3	11.2157	26.5876
4	8.5268	-3.7218
5	30.3534	5.0864
6	29.4938	-0.6784

Lets see the projected Sammons 2D onto a plane with n_cell set to 300 in second iteration.

ggplot(centroid_coordinates2, aes(x_coord, y_coord)) +
  geom_point(color = "blue") +
  labs(x = "X", y = "Y")

Figure 4: Sammons 2D Plot for 300 cells

6.2.3 Iteration 3:

lets view the projected 2D coordinates after performing sammon’s projection on the compressed data for the third iteration where we set n_cells parameter as 900. For the sake of brevity we are displaying first six rows.

hvt_torus_coordinates <-hvt.torus3[[2]][[1]][["1"]]
centroids <<- list()
  coordinates_value <- lapply(1:length(hvt_torus_coordinates), function(x){
    centroids <-hvt_torus_coordinates[[x]]
    coordinates <- centroids$pt
  })
centroid_coordinates<<- do.call(rbind.data.frame, coordinates_value)  
colnames(centroid_coordinates) <- c("x_coord","y_coord")
centroid_coordinates$Row.No <- as.numeric(row.names(centroid_coordinates)) 
centroid_coordinates <- centroid_coordinates %>% dplyr::select(Row.No,x_coord,y_coord)
centroid_coordinates3 <- centroid_coordinates %>% data.frame() %>% round(4)
Table(head(centroid_coordinates3), scroll = T, limit = 20)

Row.No	x_coord	y_coord
1	19.2964	-18.4704
2	-5.9543	10.4406
3	25.5603	0.6926
4	1.5064	-9.0975
5	18.3666	-24.9166
6	17.3898	-22.7207

Lets see the projected Sammons 2D onto a plane with n_cell set to 900 in third iteration.

ggplot(centroid_coordinates3, aes(x_coord, y_coord)) +
  geom_point(color = "blue") +
  labs(x = "X", y = "Y")

Figure 5: Sammons 2D Plot for 900 cells

6.3 Step 3: Tessellation

The deldir package computes the Delaunay triangulation (and hence the Dirichlet or Voronoi tessellation) of a planar point set according to the second (iterative) algorithm of Lee and Schacter. For subsequent levels, transformation is performed on the 2D coordinates to get all the points within its parent tile. Tessellations are plotted using these transformed points as centroids.plotHVT is the main function to plot hierarchical voronoi tessellation.

Now let’s try to understand plotHVT function. The parameters have been explained in detail below:

plotHVT(hvt.results, line.width, color.vec, pch1 = 21, centroid.size = 3, title = NULL, maxDepth = 1)

• hvt.results - A list containing the output of the HVT function which has the details of the tessellations to be plotted.

• line.width - A vector indicating the line widths of the tessellation boundaries for each layer.

• color.vec - A vector indicating the colors of the tessellations boundaries at each layer.

• pch1 - Symbol type of the centroids of the tessellations (parent levels). Refer points (default = 21).

• centroid.size - Size of centroids of first level tessellations (default = 3).

• title - Set a title for the plot (default = NULL).

• maxDepth - An integer indicating the number of levels. (default = NULL)

6.3.1 Iteration 1:

To enhance visualization, let’s generate a plot of the Voronoi tessellation for the first iteration where we set n_cells parameter as 100. This plot will provide a visual representation of the Voronoi regions corresponding to the data points, aiding in the analysis and understanding of the data distribution.

muHVT::plotHVT(
  hvt.torus,
  line.width = c(0.4),
  color.vec = c("#141B41"),
  centroid.size = 0.6,
  maxDepth = 1
)

Figure 6: The Voronoi tessellation for layer 1 shown for the 100 cells in the dataset ’torus’

6.3.2 Iteration 2:

Now, let’s plot the Voronoi tessellation for the second iteration where we set n_cells parameter to 300.

muHVT::plotHVT(
  hvt.torus2,
  line.width = c(0.4),
  color.vec = c("#141B41"),
  centroid.size = 0.6,
  maxDepth = 1
)

Figure 7: The Voronoi tessellation for layer 1 shown for the 300 cells in the dataset ’torus’

6.3.3 Iteration 3:

Now, let’s plot the Voronoi tessellation again, for the third iteration where we set n_cells parameter to 900.

muHVT::plotHVT(
  hvt.torus3,
  line.width = c(0.4),
  color.vec = c("#141B41"),
  centroid.size = 0.6,
  maxDepth = 1
)

Figure 8: The Voronoi tessellation for layer 1 shown for the 900 cells in the dataset ’torus’

From the presented plot, the inherent structure of the donut can be easily observed in the two-dimensional space

We will now overlay all the features as heatmap over the Voronoi Tessellation plot for better visualization and identification of patterns, trends, and variations in the data.

Heat Maps

Let’s have look at the hvtHmap function which we will use to overlay a variable as heatmap.

hvtHmap(hvt.results, dataset, child.level, hmap.cols, color.vec ,line.width, palette.color = 6)

• hvt.results - A list of results obtained from the HVT function.

• dataset - A dataframe containing the variables to overlay as a heatmap. The user can pass an external dataset or the dataset that was used to perform hierarchical vector quantization. The dataset should have the same number of points as the dataset used to perform hierarchical Vector Quantization in the HVT function.

• child.level - A number indicating the level for which the heat map is to be plotted.

• hmap.cols - The column number of column name from the dataset indicating the variables for which the heat map is to be plotted. To plot the quantization error as heatmap, pass 'quant_error'. Similarly to plot the no of points in each cell as heatmap, pass 'no_of_points' as a parameter.

• color.vec - A color vector such that length(color.vec) = child.level (default = NULL).

• line.width - A line width vector such that length(line.width) = child.level (default = NULL).

• palette.color - A number indicating the heat map color palette. 1 - rainbow, 2 - heat.colors, 3 - terrain.colors, 4 - topo.colors, 5 - cm.colors, 6 - BlCyGrYlRd (Blue,Cyan,Green,Yellow,Red) color (default = 6).

• show.points - A boolean indicating whether the centroids should be plotted on the tessellations (default = FALSE).

Now let’s plot the Voronoi Tessellation with the heatmap overlaid for all the features in the torus data for better visualization and interpretation of data patterns and distributions.

The heatmaps displayed below provides a visual representation of the spatial characteristics of the torus, allowing us to observe patterns and trends in the distribution of each of the features (n,X,Y and Z). The sheer green shades highlight regions with higher coordinate values in each of the heatmaps, while the indigo shades indicate areas with the lowest coordinate values in each of the heatmaps. By analyzing these heatmaps, we can gain insights into the variations and relationships between each of these features within the torus structure.

muHVT::hvtHmap(
  hvt.torus3,
  torus_df,
  child.level = 1,
  hmap.cols = "n",
  line.width = c(0.4),
  color.vec = c("#141B41"),
  palette.color = 6,
  centroid.size = 0.8,
  show.points = T,
  quant.error.hmap = 0.1,
  n_cells.hmap = 15
)

Figure 9: The Voronoi tessellation for layer 1 and number of cells 900 with the heat map overlaid for No. of entities in each cell in the ’torus’ dataset

muHVT::hvtHmap(
  hvt.torus3,
  torus_df,
  child.level = 1,
  hmap.cols = "x",
  line.width = c(0.4),
  color.vec = c("#141B41"),
  palette.color = 6,
  centroid.size = 0.8,
  show.points = T,
  quant.error.hmap = 0.1,
  n_cells.hmap = 15
)

Figure 10: The Voronoi tessellation for layer 1 and number of cells 900 with the heat map overlaid for variable x in the ’torus’ dataset

muHVT::hvtHmap(
  hvt.torus3,
  torus_df,
  child.level = 1,
  hmap.cols = "y",
  line.width = c(0.4),
  color.vec = c("#141B41"),
  palette.color = 6,
  centroid.size = 0.8,
  show.points = T,
  quant.error.hmap = 0.1,
  n_cells.hmap = 15
)

Figure 11: The Voronoi tessellation for layer 1 and number of cells 900 with the heat map overlaid for variable y in the ’torus’ dataset

muHVT::hvtHmap(
  hvt.torus3,
  torus_df,
  child.level = 1,
  hmap.cols = "z",
  line.width = c(0.4),
  color.vec = c("#141B41"),
  palette.color = 6,
  centroid.size = 0.8,
  show.points = T,
  quant.error.hmap = 0.1,
  n_cells.hmap = 15
)

Figure 12: The Voronoi tessellation for layer 1 and number of cells 900 with the heat map overlaid for variable z in the ’torus’ dataset

6.4 Step 4: Prediction(predictHVT)

Raw Testing Dataset

Lets have a look at our randomly selected test dataset containing (3000 points) before we pass it to predictHVT function for scoring.

Table(head(test_dataset1))

Row.No	x	y	z
6	-1.4824	0.9229	0.9672
10	0.7920	-1.3482	-0.8998
12	-2.3787	1.7986	-0.1878
17	-0.8428	-0.5436	0.0755
20	-2.6487	-0.5745	0.7040
23	-1.1130	-0.6516	-0.7040

However, let’s try to comprehend the predictHVT function first before moving on

predictHVT(data,
                  hvt.results,
                  hmap.cols = NULL,
                  child.level = 1,
                  ...)

The important parameters for the function predictHVT are as below:

• data - A dataframe containing the test dataset. The dataframe should have all the variable(features) used for training. The variables from this dataset can also be used to overlay as heatmap.

• hvt.results - A list of hvt.results obtained from the HVT function while performing hierarchical vector quantization on training data. The list containes detailed information about the hierarchical vector quantized data along with a summary section containing no of points, Quantization Error and the centroids for each cell, as per the n_cells given to the HVT() function.

• hmap.cols - The column number of column name from the dataset indicating the variables for which the heat map is to be plotted. A heatmap won’t be plotted if NULL is passed (Default = NULL)

• child.level - A number indicating the level for which the heat map is to be plotted (Only used if hmap.cols is not NULL) Each level represents a different level of clustering or partitioning of the data.

• normalize - A logical value indicating if the columns in your dataset should be normalized. Basically it is a technique that scales the values of each variable to have a mean of 0 and a standard deviation of 1.. Default value is TRUE.

• distance_metric - It specifies the type of distance measurement used to calculate similarity or dissimilarity between data points. It can be set to “Euclidean” (default) for straight-line distance or “Manhattan” for the sum of absolute differences between coordinates.

• error_metric - It specifies the error metrics to be used for evaluating the performance of the model. It can be “mean” or “max”. mean is selected by default.

• yVar - Name of the dependent variable(s)

• ... - color.vec and line.width can be passed from here

Now once we have built the model, let us try to predict using our test dataset which cell and which level each point belongs to.

set.seed(240)
predictions_torus <- muHVT::predictHVT(
  testTorus,
  hvt.torus3,
  child.level = 1,
  line.width = c(1.2),
  color.vec = c("#141B41"),
  quant.error.hmap = 0.1,
  n_cells.hmap = 900,
  normalize = F
)

Let’s see which cell and level each point belongs to and check the mean absolute difference. For the sake of brevity, we will only show the first 10 rows

data1 <- test_dataset
data1$Row.No <- row.names(test_dataset)
data1 <- data1 %>% dplyr::select(Row.No,x,y,z)
rownames(data1) <- NULL
colnames(data1) <- c("Row.No","x_act","y_act","z_act")
data2 <- predictions_torus[["scoredPredictedData"]]
data2 <- data2 %>% dplyr::select(Cell.ID,x,y,z)
colnames(data2) <- c("Cell.ID","x_pred","y_pred","z_pred")
combined <- cbind(data1,data2)
combined$diff <- rowMeans(abs(combined[, c("x_act", "y_act", "z_act")] - combined[, c("x_pred", "y_pred", "z_pred")]))
options(scipen = 999)
combined %>% head(100) %>% 
  as.data.frame() %>%
  Table(scroll = T, limit = 10)

Row.No	x_act	y_act	z_act	Cell.ID	x_pred	y_pred	z_pred	diff
6	-1.4823709	0.9228529	0.9672467	723	-1.4824	0.9229	0.9672	0.0000410
10	0.7920450	-1.3482111	-0.8997781	252	0.7920	-1.3482	-0.8998	0.0000260
12	-2.3787465	1.7986402	-0.1878163	900	-2.3787	1.7986	-0.1878	0.0000344
17	-0.8427718	-0.5435588	0.0755262	558	-0.8428	-0.5436	0.0755	0.0000318
20	-2.6486525	-0.5744624	0.7039659	837	-2.6487	-0.5745	0.7040	0.0000397
23	-1.1130327	-0.6516259	-0.7039507	628	-1.1130	-0.6516	-0.7040	0.0000360
28	0.7520403	-2.6043863	0.7034024	140	0.7520	-2.6044	0.7034	0.0000188
30	-1.6755282	2.3358857	0.4847097	859	-1.6755	2.3359	0.4847	0.0000174
33	-1.6466922	0.4011827	-0.9523068	719	-1.6467	0.4012	-0.9523	0.0000106
34	0.7930278	2.4427927	0.8228262	458	0.7930	2.4428	0.8228	0.0000204

hist(combined$diff, breaks = 20, col = "blue", main = "Mean Absolute Difference", xlab = "Difference")

Figure 13: Mean Absolute Difference

7 Example II: muHVT with the Personal Computer dataset

Data Understanding

In this section, we will use the Prices of Personal Computers dataset. This dataset contains 6259 observations and 10 features. The dataset observes the price from 1993 to 1995 of 486 personal computers in the US. The variables are price, speed, ram, screen, cd, etc. The dataset can be downloaded from here.

In this example, we will compress this dataset by using hierarchical VQ via k-means and visualize the Voronoi Tessellation plots using Sammons projection. Later on, we will overlay all the variables as a heatmap to generate further insights.

Here, we load the data and store into a variable computers.

set.seed(240)
# Load data from csv files
computers <- read.csv("https://raw.githubusercontent.com/Mu-Sigma/muHVT/master/vignettes/sample_dataset/Computers.csv")

Raw Personal Computers Dataset

The Computers dataset includes the following columns:

• Price: The price of 486 PCs in US dollars.
• Speed: The clock speed of the PCs in MHz.
• HD: The size of the hard drive in MB.
• RAM: The size of RAM in MB.
• Screen: The size of the screen in inches.
• Ads: The number of listings for 486 PCs’ prices for each month.

Let’s explore the Personal Computers Dataset containing (6259 points). For the sake of brevity we are displaying first six rows.

# Quick peek
Table(head(computers), scroll = T, limit = 20)

Now, let us check the structure of the data and analyse its summary.

str(computers)
#> 'data.frame':    6259 obs. of  11 variables:
#>  $ X      : int  1 2 3 4 5 6 7 8 9 10 ...
#>  $ price  : int  1499 1795 1595 1849 3295 3695 1720 1995 2225 2575 ...
#>  $ speed  : int  25 33 25 25 33 66 25 50 50 50 ...
#>  $ hd     : int  80 85 170 170 340 340 170 85 210 210 ...
#>  $ ram    : int  4 2 4 8 16 16 4 2 8 4 ...
#>  $ screen : int  14 14 15 14 14 14 14 14 14 15 ...
#>  $ cd     : chr  "no" "no" "no" "no" ...
#>  $ multi  : chr  "no" "no" "no" "no" ...
#>  $ premium: chr  "yes" "yes" "yes" "no" ...
#>  $ ads    : int  94 94 94 94 94 94 94 94 94 94 ...
#>  $ trend  : int  1 1 1 1 1 1 1 1 1 1 ...

summary(computers)
#>        X            price          speed              hd        
#>  Min.   :   1   Min.   : 949   Min.   : 25.00   Min.   :  80.0  
#>  1st Qu.:1566   1st Qu.:1794   1st Qu.: 33.00   1st Qu.: 214.0  
#>  Median :3130   Median :2144   Median : 50.00   Median : 340.0  
#>  Mean   :3130   Mean   :2220   Mean   : 52.01   Mean   : 416.6  
#>  3rd Qu.:4694   3rd Qu.:2595   3rd Qu.: 66.00   3rd Qu.: 528.0  
#>  Max.   :6259   Max.   :5399   Max.   :100.00   Max.   :2100.0  
#>       ram             screen           cd               multi          
#>  Min.   : 2.000   Min.   :14.00   Length:6259        Length:6259       
#>  1st Qu.: 4.000   1st Qu.:14.00   Class :character   Class :character  
#>  Median : 8.000   Median :14.00   Mode  :character   Mode  :character  
#>  Mean   : 8.287   Mean   :14.61                                        
#>  3rd Qu.: 8.000   3rd Qu.:15.00                                        
#>  Max.   :32.000   Max.   :17.00                                        
#>    premium               ads            trend      
#>  Length:6259        Min.   : 39.0   Min.   : 1.00  
#>  Class :character   1st Qu.:162.5   1st Qu.:10.00  
#>  Mode  :character   Median :246.0   Median :16.00  
#>                     Mean   :221.3   Mean   :15.93  
#>                     3rd Qu.:275.0   3rd Qu.:21.50  
#>                     Max.   :339.0   Max.   :35.00

Let us first split the data into train and test. We will randomly select 80% of the data as train and remaining as test.

num_rows <- nrow(computers)
set.seed(123)
train_indices <- sample(1:num_rows, 0.8 * num_rows)
trainComputers <- computers[train_indices, ]
testComputers <- computers[-train_indices, ]

K-means is not suitable for factor variables as the sample space for factor variables is discrete. A Euclidean distance function on such a space isn’t really meaningful. Hence, we will delete the factor variables(X, cd, multi, premium, trend) in our dataset.

trainComputers <-
  trainComputers %>% dplyr::select(-c(X, cd, multi, premium, trend))
testComputers <-
  testComputers %>% dplyr::select(-c(X, cd, multi, premium, trend))

Raw Training Dataset

Now, lets have a look at the randomly selected raw training dataset containing (5007 data points). For the sake of brevity we are displaying first six rows.

trainComputers_data <- trainComputers %>% as.data.frame() %>% round(4)
trainComputers_data$Row.No <- as.numeric(row.names(trainComputers_data))
trainComputers_data <- trainComputers_data %>% dplyr::select(Row.No,price,speed,hd,ram,screen,ads)
row.names(trainComputers_data) <- NULL
Table(head(trainComputers_data))

Raw Testing Dataset

Now, lets have a look at the randomly selected raw testing dataset containing (1252 data points). For the sake of brevity we are displaying first six rows.

#testComputers <- scale(testComputers, center = scale_attr$`scaled:center`, scale = scale_attr$`scaled:scale`) 
testComputers_data <- testComputers %>% as.data.frame() %>% round(4)
testComputers_data$Row.No <- as.numeric(row.names(testComputers_data))
testComputers_data <- testComputers_data %>% dplyr::select(Row.No,price,speed,hd,ram,screen,ads)
rownames(testComputers_data) <- NULL
Table(head(testComputers_data))

As we are familiar with the structure of the computers data, we will now follow the following steps to get the predictions using the Computers dataset.

set.seed(240)
hvt.results <- list()
hvt.results <- muHVT::HVT(trainComputers,   
                          n_cells = 440,
                          depth = 1,
                          quant.err = 0.2,
                          projection.scale = 10,
                          normalize = T,
                          distance_metric = "L1_Norm",
                          error_metric = "max",
                          quant_method = "kmeans",
                          diagnose = F)

compressionSummaryTable(hvt.results[[3]]$compression_summary)

summaryTable(hvt.results[[3]]$summary)